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ihe  goal  of  the  program  on  Applications  of  Molecular  Dynemics  to  the 
Study  of  Shock-Induced  Chemistry  is  to  develop  and  apply  molecular  dynamics 
techniques  to  problems  in  shock-induced  chemistry  in  reactive  and  non¬ 
reactive  materials.  The  three  major  thrusts  of  the  work  proposed  are:  (1) 
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To  develop  phenomenological  submodels  of  molecular  Interactions  that  can  b 
used  In  full  molecular  dynamics  calculations.  In  particular,  to  develop 
models  which  include  known  quantum  mechanical  properties  of  atoms  and 
molecules,  (2)  to  develop  a  model  for  doing  many  particle  Interactions  or 
vector  computer,  and  (3)  to  apply  the  resulting  numerical  models  to  probl j 
of  shocks  in  condensed  phase  materials. 

The  first  part  of  this  research  has  been  done  in  close  collaboration 
with  Dr.  Hersch  Rabitz  at  Princeton  University  and  Dr.  Robert  Wyatt  at  th> 
University  of  Texas  in  Austin.  Recently,  Dr.  Boyd  Waite  from  the  U.S.  Na ; 
Academy  has  also  been  participating.  The  work  has  been  coordinated  with 
Special  Focus  Programs  in  Energetic  Materials  at  NRL  and  ONR.  The  work  a 
NRL  is  partially  supported  by  a  6.1  NRL  research  program  in  Computational 
Physics. 

To  date,  we  have  concentrated  on  the  first  two  aspects  of  the  resear 
We  have  worked  with  Drs.  Wyatt  and  Rabitz  in  deciding  how  to  set  up  a  tes 
problem  that  could  be  used  to  test  the  way  in  which  quantum  mechanical 
effects  can  be  phenomenologically  modelled  in  a  classical  representation. 
Also,  we  have  begun  development  and  testing  of  a  new  vectorized  algorithm 
for  performing  molecular  dynamics  simulations  that  will  take  advantage  of 
the  abilities  of  a  CRAY.  The  progress  on  these  areas  is  summarized  in  th 
report. 


Accession 
NT  IS 
DTIC  ?' 

Un:,n\  .. 

Just  !  .  ■! 


L 1  :■  I  r ' ,  •  -  1  v  ■.  / 

Av.;il  .1  ti  :  l  '.'.y  Codes  / 
;-..U  -.../or  t 

)ist  >  >i‘t.  lt*l 


1984  ANNUAL  REPORT  FOR  THE  PROJECT 


APPLICATIONS  OF  MOLECULAR  DYNAMICS  TO  THE  STUDY  OF 
SHOCK-INDUCED  CHEMISTRY 


Principal  Investigator:  Dr.  Elaine  S.  Oran 

Laboratory  for  Computational  Physics 
Code  4040 

Washington,  D.C.  20375 

Dr.  Elaine  S.  Oran 
Dr.  Jay  P.  Boris 
Dr.  Michael  Page 

Dr.  Sam  Lambrakos  (NRC-Postdoctoral  Research  Associate) 

Dr.  Boyd  Waite  (U.S.  Naval  Academy) 

Ms.  Lauree  Shampine  (summer  student,  junior  at  Bryn  Mawr  College) 
Ms.  Liz  Gold  (summer  student,  freshman  at  Princeton  University) 

The  goal  of  this  program  is  to  develop  and  apply  molecular  dynamics 
techniques  to  problems  in  shock-induced  chemistry  in  reaotive  and  non¬ 
reactive  materials.  The  three  major  thrusts  of  the  work  proposed  are: 

1.  To  develop  phenomenological  submodels  of  molecular  interactions 
that  can  be  used  in  full  molecular  dynamics  calculations.  In  particular, 
to  develop  models  which  include  known  quantum  mechanical  properties  of 
atoms  and  molecules. 

2.  To  develop  a  model  for  doing  many  particle  interactions  on  a 
vector  computer, 

3.  To  apply  the  resulting  numerical  models  to  problems  of  shooks  in 
condensed  phase  materials. 

The  first  part  of  this  research  has  been  done  in  close  collaboration 
with  Dr.  Hersch  Rabitz  at  Princeton  University  and  Dr.  Robert  Wyatt  at  the 
University  of  Texas  in  Austin,  Recently,  Dr.  Boyd  Waite  from  the  U.S.  Naval 
Academy  has  also  been  participating.  The  work  has  been  coordinated  with  the 
Special  Focus  Programs  in  Energetic  Materials  at  NRL  and  ONR.  The  work  at 
NRL  is  partially  supported  by  a  6.1  NRL  research  program  in  Computational 
Physics 
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CURRENT  STATUS  OF  RESEARCH 


To  date,  we  have  concentrated  on  the  first  two  aspects  of  the  research. 

We  have  worked  with  Drs.  Wyatt  and  Rabitz  in  deciding  how  to  set  up  a  test 
problem  that  could  be  used  to  test  the  way  in  which  quantum  mechanical  effects 
can  be  phenomenologically  modelled  in  a  classical  representation.  Also,  we 
have  begun  development  and  testing  of  a  new  vectorized  algorithm  for  performing 
molecular  dynamics  simulations  that  will  take  advantage  of  the  abilities  of  a 
CRAY. 

Development  of  Submodels  for  Molecular  Dynamics 

Our  goal  is  to  create  models  of  molecular  Interactions  that  represent 
some  of  the  known  quantum  mechanical  effects  that  do  not  appear  in  a  purely 
classical  representation.  Such  models  could  then  be  used  in  molecular  dynamics 
calculations  and  would  be  one  way  of  bridging  the  gap  between  classical  and 
quantum  calculations.  Selected  teat  problems  are  being  solved  classically, 
semi-classically ,  and  quantum  mechanically  in  order  to  develop  these  models. 

After  extensive  discussions  we  decided  on  a  test  problem  in  which  there 
are  two  collinear  quantum  harmonic  osolllators,  l.e.,  four  particles,  for 
which  we  nominally  chose  the  masses  and  force  constants  to  correspond  to  the 
oxygen  and  nitrogen  molecules.  The  closest  oxygen  and  nitrogen  atom  Interact 
through  a  Leonard- Jones  potential.  This  problem  was  solved  three  ways: 
classically  (NRL),  semi-classically  (NRL),  and  quantum  mechanically  (Texas). 

The  total  energy  of  the  system  was  to  be  large  enough  so  that  energy  oould  go 
into  the  harmonic  oscillator  states  and  there  would  still  be  enough  energy 
left  for  an  unbound  Leonard- Jones  collision.  Below  we  describe  the  current 
status  of  each  of  these  calculations. 
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(a)  Quantum  Mechanical  Calculations 

A  code  to  do  the  test  problen  quantum  mechanically  was  set  up  at 
the  University  of  Texas  by  Dr.  Wyatt  and  his  student,  Mr.  Don  Miller. 
From  this  code  ve  have  extracted  those  portions  required  for  the  semi- 
classical  calculation,  i.e.  the  potential  energy  matrix.  They  have  run 


the  code  for  the  test  problem. 

As  the  code  is  now  written,  direct  comparisons  can  be  made  of  the 
final  state  of  the  two  molecules,  given  an  initial  state.  However,  we 


are  Interested  in  comparing  the  results  along  the  course  of  a  collision 
in  order  to  develop  models  to  update  the  states  of  many  particles  at  fixed 
short  Intervals  in  the  molecular  dynamics  simulations.  We  have  not  yet 
found  a  way  to  do  this  with  their  current  formulation  without  a  great  deal 
of  effort  on  their  part.  Thus  for  now,  we  can  only  compare  their  final 
answers  to  the  semi-classical  and  classical  calculations. 

(b)  Semi-Classical  Calculations 

There  are  three  aspects  of  the  work  on  the  semi-classical  formulation 
that  we  describe  here:  (1)  finding  the  best  way  to  solve  the  ordinary 
differential  equations,  (2)  how  to  make  the  solutions  conserve  energy,  i.e., 
how  to  handle  the  quantum-classical  feedback  problem;  and  (3)  the  recent 
Interest  of  Dr,  Boyd  Waite  from  Annapolis  in  solving  the  test  problem  by  the 
methods  pioneered  by  Dr.  William  Miller  from  the  University  of  California, 
Berkeley. 

(1)  A  number  of  integration  methods  for  the  ordinary  differential 
equations  were  tried.  It  was  quickly  determined  the  even  a  sophisticated 
fourth  order  Runge-Kutta  was  not  good  enough  at  the  turn-around  point  of  a 
collision.  We  were  able  to  solve  this  problem  by  using  a  trick  which 
we  had  already  used  very  successfully  in  the  classical  calculations  described 
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below:  use  a  leap-frog  method,  which  la  second  order,  for  the  classical 
variables.  The  difference  here  was  that  we  still  used  the  Runge-Kutta  for 
the  quantum  mechanical  variables.  This  gives  excellent  results,  due  to 
the  reversible  property  of  a  leap-frog  method.  It  also  allows  us  to  take 
larger  timesteps  and  do  the  calculation  much  more  quickly.  The  variable 
timestep  Runge-Kutta  required  an  Inordinate  amount  of  time,  especially  at 
the  turn-around  points  and  lack  of  microscopic  reversibility  led  to  a 
nonphysical  drift  in  the  total  energy  of  the  system. 

(2)  In  the  first  semi-classical  calculations  we  performed,  no 
feedback  was  allowed  to  the  classical  trajectory  from  anything  happening 
to  the  quantum  mechanical  variables.  The  molecules  approached  each  other, 
and  the  kinetic  energy  of  the  collision  was  converted  into  internal  energy 
of  the  molecules  through  the  Lennard-Jones  potential.  However,  no  mechanism 
existed  in  the  formulation  to  slow  down  the  collision  as  a  result  of  energy 
being  extracted.  When  the  molecules  separated,  no  mechanism  existed  to 

take  energy  from  the  internal  states  of  the  molecules  and  put  it  into  kinetic 
energy  of  the  system.  Thus  the  result  showed  non-conservation  of  energy. 

Our  solution  was  to  devise  a  feedback  mechanism  which  introduced  this  vibra¬ 
tion-translation  energy  transfer  in  such  a  way  that  the  energy  conservation 
constraint  remained  satisfied.  The  idea  is  to  monitor  the  expectation  value 
of  the  unperturbed  harmonic  oscillator  Hamiltonians  and  adjust  the  classical 
trajectory  accordingly.  This  procedure  gives  a  classical  path  whioh  is  a 
weighted  average  of  the  individual  atate-to-state  trajectories  at  each  point. 
These  calculations  are  complete  now,  and  must  be  compared  with  the  the  quantum 
mechanical  and  classical  solutions. 

(3)  Dr.  Boyd  Halts  from  the  Chemistry  Department  at  Annapolis  has  done 
semi-classical  ealoulatlons  using  the  methods  developed  by  Dr.  William  Miller 


from  the  University  of  California,  Berkeley.  Dr.  Waite  is  now  applying  his 
methods  to  the  test  problem  we  are  doing. 

(c)  Classical  Calculations 

A  rather  general  classical  code  was  written  that  will  treat  an 
arbitrary  number  of  particles,  each  pair  of  which  will  interact  with  an 
arbitrarily  specified  potential.  This  is  not  meant  to  be  a  full  molecular 
dynamics  model,  but  something  which  can  be  used  as  a  test-bed  for  various 
potentials  and  models  to  be  developed  for  the  large  scale  molecular 
dynamics  model  described  below.  This  model  was  configured  for  the 
four-particle  test  problem  and  a  number  of  preliminary  tests  have  been  run. 

The  first  tests  involved  integration  method  tests.  Here  we 
thoroughly  tested  the  leap-frog  method  and  found  that  it  was  extremely 
accurate  and  did  an  excellent  job  at  conserving  energy.  The  method 
was  best  when  we  kept  the  timestep  constant,  which  is  required  to 
maintain  the  reversibility  property.  We  have  also  devised  an  algorithm 
which  keeps  the  distance  between  two  atoms  in  a  molecule  fixed  and  still 
retains  the  necessary  microsoopio  reversibility  in  the  integration. 

Calculations  were  then  done  to  test  the  validity  of  the  expansion 
of  the  Lennard-Jones  potential  developed  at  the  University  of  Texas  by 
comparing  a  range  of  classical  calculations  using  it  and  the  full 
potential.  The  expansion  is  excellent  and  for  the  worst  cases  tested, 
differences  appeared  only  in  the  third  or  fourth  decimal  place. 

We  then  performed  the  classical  calculations  of  the  test  problem 
to  oompare  with  the  seml-olasslcal  and  quantum  mechanical  results. 

This  involves  many  thousands  of  oolllslon  calculations  with  this  code, 
each  with  varying  initial  energies  and  phases.  There  are  a  number  of 
different  ways  to  perform  the  statistics  on  these  calculations.  The 


different  ways  involve  different  procedures  for  choosing  initial  energies 
and  phases  and  binning  the  results.  Dr.  Waite  consulted  with  us  on  how 
to  fornulate  the  statistics  on  the  classical  problen.  We  have  also  tested 
our  classical  calculation  code  against  his  as  a  benchmark,  and  currently 
his  programs  are  also  being  used  to  solve  the  problem  classically. 

Molecular  Dynamics  Calculations 

We  are  currently  developing  a  three-dimensional  vectorized  molecular 
dynamics  model  that  will  host  the  submodels  described  above.  Performing 
accurate  enough  calculations  with  enough  particles  to  get  reasonable 
statistics  is  the  heart  of  the  molecular  dynamics  problem.  One  of  the 
key  problems  here  is  computing  the  interaction  of  each  particle  with  all  of 
the  particles  around  it.  Among  N  particles,  there  are  N-squared  interactions 
to  be  considered.  The  computer  costs  to  perform  all  of  these  calculations 
is  exhorbitant  without  a  suitable  trick. 

This  problem  has  led  us  to  devise  a  three-dimensional  nearest  neighbor 
algorithm  whose  cost  scales  as  the  number  M  of  independant  objects,  not  the 
square  of  N  as  with  previous  algorithms.  The  algorithm  is  based  on  a 
dynamic  monotonic  logical  grid  whioh  is  constructed  so  that  nearby  data  in 
the  logical  grid  is  automatically  guaranteed  to  correspond  to  points  which 
are  nearby  in  real  space.  It  ezeoutes  in  small  array  processors  and  also 
partitions  directly  to  take  advantage  of  multitasking  asynchronous 
parallelism  in  distributed  processing  systems.  We  have  shown  that  the 
present  test  model  executes  extremely  fast  on  a  problem  Involving  512 
particles,  and  we  believe  that  this  is  an  order  of  magnitude  faster  than 
any  other  current  existing  method.  The  relative  gains  will  be  even  greater 
for  larger  systems  because  of  the  favorable  soallng. 


PROPOSED  RESEARCH 


In  FI’ 85  we  will 

1.  Complete  the  comparisons  of  the  solution  of  the  test  problem 
which  was  solved  classically,  semi-classlcally,  and  quantum  mechanically. 

2.  Begin  development  of  phenomenological  models  for  including  quantum 
mechanical  effects  in  classical  molecular  dynamics  calculations. 

3.  Set  up  the  molecular  dynamics  calculation  with  the  new  monotonic 
logical  grid  algorithm.  Test  it  on  a  standard  condensation  problem  after 
a  number  of  suitable  tests  of  the  new  griddlng  algorithms.  Include  a 
framework  for  Inserting  phenomenological  quantum  models. 

PUBLICATIONS  and  PRESENTATIONS 

A  Comparison  of  Quantal,  Classical,  and  Semidasslcal  Descriptions  of 
a  Model  Collinear  Inelastic  Diatomic  Collision,  currently  being  completed, 
M.Page,  D.  Miller,  R.  Vyatt,  H.  Rabitz,  E.  Oran,  and  B.  Vaite. 

A  Vectorized  "Nearest-Neighbor"  Algorithm  of  Order  N  Using  a  Monotonic 
Logical  Grid,  currently  being  completed,  J.  Boris. 

A  Vectorized  "Nearest-Neighbor"  Algorithm  of  Order  N  Using  a  Monotonic 
Logical  Grid,  J.  Boris,  Meeting  of  the  American  Physical  Society, 

Boston,  October,  1934. 
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